home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / deriv.pro < prev    next >
Text File  |  1997-07-08  |  2KB  |  68 lines
  1. ; $Id: deriv.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1984-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;
  6.  
  7. Function Deriv, X, Y
  8. ;+
  9. ; NAME:
  10. ;    DERIV
  11. ;
  12. ; PURPOSE:
  13. ;    Perform numerical differentiation using 3-point, Lagrangian 
  14. ;    interpolation.
  15. ;
  16. ; CATEGORY:
  17. ;    Numerical analysis.
  18. ;
  19. ; CALLING SEQUENCE:
  20. ;    Dy = Deriv(Y)         ;Dy(i)/di, point spacing = 1.
  21. ;    Dy = Deriv(X, Y)    ;Dy/Dx, unequal point spacing.
  22. ;
  23. ; INPUTS:
  24. ;    Y:  Variable to be differentiated.
  25. ;    X:  Variable to differentiate with respect to.  If omitted, unit 
  26. ;        spacing for Y (i.e., X(i) = i) is assumed.
  27. ;
  28. ; OPTIONAL INPUT PARAMETERS:
  29. ;    As above.
  30. ;
  31. ; OUTPUTS:
  32. ;    Returns the derivative.
  33. ;
  34. ; COMMON BLOCKS:
  35. ;    None.
  36. ;
  37. ; SIDE EFFECTS:
  38. ;    None.
  39. ;
  40. ; RESTRICTIONS:
  41. ;    None.
  42. ;
  43. ; PROCEDURE:
  44. ;    See Hildebrand, Introduction to Numerical Analysis, Mc Graw
  45. ;    Hill, 1956.  Page 82.
  46. ;
  47. ; MODIFICATION HISTORY:
  48. ;    Written, DMS, Aug, 1984.
  49. ;-
  50. ;
  51.  
  52. on_error,2              ;Return to caller if an error occurs
  53. n = n_elements(x)
  54. if n lt 3 then message, 'Parameters must have at least 3 points'
  55.  
  56. if n_params(0) ge 2 then begin
  57.     if n ne n_elements(y) then message,'Vectors must have same size'
  58.     d = float(shift(y,-1) - shift(y,1))/(shift(x,-1) - shift(x,1))
  59.     d[0] = (-3.0*y[0] + 4.0*y[1] - y[2])/(x[2]-x[0])
  60.     d[n-1] = (3.*y[n-1] - 4.*y[n-2] + y[n-3])/(x[n-1]-x[n-3])
  61.    end else begin
  62.     d = (shift(x,-1) - shift(x,1))/2.
  63.     d[0] = (-3.0*x[0] + 4.0*x[1] - x[2])/2.
  64.     d[n-1] = (3.*x[n-1] - 4.*x[n-2] + x[n-3])/2.
  65.    endelse
  66. return, d
  67. end
  68.